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Abstract 

In this work we consider optical flow on evolving Riemannian 2-manifolds 
which can be parametrised from the 2-sphere. Our main motivation is to 
estimate cell motion in time-lapse volumetric microscopy images depict¬ 
ing fluorescently labelled cells of a live zebrafish embryo. We exploit the 
fact that the recorded cells float on the surface of the embryo and allow 
for the extraction of an image sequence together with a sphere-like sur¬ 
face. We solve the resulting variational problem by means of a Galerkin 
method based on vector spherical harmonics and present numerical results 
computed from the aforementioned microscopy data. 


1 Introduction 

Motion estimation is a fundamental problem in image analysis and computer 
vision. An important task within is optical flow computation. It is concerned 
with the inference of a vector field describing the displacements of brightness 
patterns, such as moving objects, in a sequence of images. Ever since the seminal 
work of Horn and Schunck [18] a variety of reliable and efficient methods have 
been proposed and successfully applied in a wide number of fields. 

Primarily, optical flow is computed in the plane. However, it is readily gen¬ 
eralised to non-Euclidean settings allowing, for instance, for cell motion analysis 
in time-lapse microscopy data. It has been only recently that high-resolution 
observations of biological model organisms such as the zebrafish became pos- 
sibie. Despite its importance for tissue and organ formation, iittie is known 
about ceii migration and proiiferation patterns during the zebrafish’s eariy em¬ 
bryonic deveiopment [DIM]- Eiuorescence microscopy nowadays aiiows to record 
time-iapse images on the scaie of singie ceiis, see e.g. 1201123 El- Increasing 
spatiai as weii as temporai resoiutions resuit in vast amounts of data, render¬ 
ing extraction of information through visuai inspection carried out by humans 
impracticabie. Automated ceii motion estimation therefore is key to iarge-scaie 
anaiysis of such data. Opticai flow computation deiivers necessary quantitative 
methods and ieads to insights into the underiying ceiiuiar mechanisms and the 
dynamic behaviour of ceiis. See, for exampie, I1I211E3E1] and the references 
therein. 

The primary bioiogicai motivation for this work is the desire to anaiyse ceii 
motion in a iiving zebrafish during eariy embryogenesis. The data at hand depict 
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Figure 1: Frames 70 (left) and 71 (right) of the volumetric zebrafish microscopy 
images recorded during early embryogenesis. The sequence contains a total 
number of 75 frames. Fluorescence response is indicated by blue colour and is 
proportional to the observed intensity. All dimensions are in micrometer (/im). 


endodermal cells expressing a green fluorescent protein. By virtue of laser¬ 
scanning microscopy, (volumetric time-lapse) 4D images of these labelled cells 
can be recorded without capturing the background. It is known that endodermal 
cells float on a so called monolayer during early embryonic development meaning 
that they do not stack on top of each other [32]. Figure depicts two frames 
of the captured sequence, containing only the upper hemisphere of the animal 
embryo. Observe the salient formation of the cells and the noise present in the 
images. More precisely, one can see the nuclei of cells forming a round surface 


in a single layer. For more details on the microscopy data we refer to Sec. 5.1 


We exploit this situation and model this layer as an evolving surface. A nat¬ 
ural candidate for a parametrisation of such a zebrafish embryo is a sphere-like 
surfaee. It is topologically diffeomorphic to the 2-sphere and most commonly 
defined as the set of points 


{p{x)x : X e 5 ^}. 


The function p : (0, c>o) can be thought of as a radial deformation of 

and will have a dependence on time in the present paper. As a consequence, 
changes in the embryo’s geometry are attributed accordingly, albeit valid only 
during early stages of its development as cells tend to cluster subsequently. The 
main intention of this work is to conceive cell motion only on this moving 2- 
dimensional manifold. As a result we are able to reduce the spatial dimension 
of the data allowing for more efficient motion estimation in microscopy data. 
Figure depicts two frames of the surface together with images obtained by 
restriction of the volumetric microscopy data in Fig. 

In this work we model the data as a time-dependent non-negative function 
/. Its value directly corresponds to the fluorescence response of the observed 
cells. For a fixed time instant t G [0,T], the domain of / is presumed to be a 
closed surface Ait C We assume that this surface can be parametrised by 
a smooth radial map from the 2-sphere. The temporal evolution of the data / 
can then be tracked by solving an optical flow problem on this moving surface 
or, more conveniently, an equivalent problem on the round sphere. 

Traditionally, the starting point for optical flow is the assumption of constant 
brightness: a point moving along a trajectory does not change its intensity over 
time. On a moving domain A4 = {Ait}t one equivalently seeks, for every time 
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Figure 2: Depicted are frames no. 70 (left) and 71 (right) of the processed 
zebrafish microscopy sequence. Top and bottom row differ by a rotation of 180 
degrees around the xs-axis. All dimensions are in micrometer (/im). 


t G [0, T], a tangent vector field v that solves a generalised optical flow equation 

rff/ + VAl/-V = 0 (1) 

at every point x G A4, where / is the image sequence living on A4. Here, for a 
fixed time m denotes the (spatial) surface gradient, dot the standard inner 
product, and f an appropriate temporal derivative. 

The optical flow problem is ill-posed meaning that equation 0 is not uniquely 
solvable. A common approach to deal with non-uniqueness is Tikhonov regu- 
larisation, which consists of computing a minimiser of 

^a(v) = r>(v, /) + an{w). 

The first term of the sum is usually the squared norm of the left-hand side 
of 0 and, in the present article, the second term will be an Sobolev norm. 
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1.1 Contributions 


The primary concern of this article is optical flow computation on evolving 2- 
dimensional Riemannian manifolds which can be parametrised from the sphere. 
Motivated by the aforementioned zebrafish microscopy data we consider closed 
surfaces for which the mapping 

{t,x) ^ p{t,x)x, X e (2) 

is a diffeomorphism between the 2-sphere and Aft for every time t G [0,T]. As 
a prototypical example we restrict ourselves to radially parametrised surfaces 
as they suit quite naturally to the given data. 

The contributions of this work are as follows. First, we give a variational 
formulation of optical flow on 2-dimensional closed Riemannian manifolds. We 
assume a dependence on time and speak of evolving surfaces. The main idea is 
to solve the problem by a Galerkin method in a finite-dimensional subspace of 
an appropriate (vectorial) Sobolev space. We take advantage of the fact that 
tangential vector spherical harmonics form a complete orthonormal system for 
1/^(5^,T5^). The sought vector field is thus uniquely determined when ex¬ 
panded in terms of the pushforward—by means of the differential of —of 
these functions. From that we arrive at a minimisation problem over where 
n is the dimension of the finite-dimensional space, and state the optimality con¬ 
ditions. They can be written purely in terms of spherical quantities and solved 
on the 2-sphere. To this end, we use a standard polyhedral approximation and 
locally interpolate spherical functions by piecewise quadratic polynomials. For 
numerical integration we employ appropriate quadrature rules on the approxi¬ 
mated sphere. 

Second, to obtain the smooth sphere-like surface, which is described by the 
map from the observed microscopy data, we formulate another variational 
problem on the sphere. The problem is essentially surface interpolation with 

Sobolev seminorm regularisation. Approximate cell centres serve as sample 
points of the surface. In particular, our microscopy data are supported only on 
the upper hemisphere, see Figs. and Scalar spherical harmonics are the 
appropriate choice for the numerical solution of the surface fitting problem, as 
they provide great flexibility with respect to the chosen space . 

Finally, we present numerical experiments on the basis of the mentioned cell 
microscopy data of a live zebrafish. To this end we compute an approximation 
of the sphere-shaped embryo and obtain a sequence of images living on this 
moving surface. Eventually, we solve for the optical flow and present the results 
in a visually adequate manner. 

1.2 Related Work 

The first variational formulation of optical flow is commonly attributed to Horn 
and Schunck m- They attempted to compute a displacement field in by min¬ 
imising a Tikhonov-regularised energy functional. It favours spatially regular 
vector fields by penalising its squared Sobolev seminorm. For introductory 
material on the subject we refer to HE] and to m for a survey on various 
optical flow functionals. Well-posedness of the aforementioned energy was first 
shown by Schnorr [35]. Moreover, there the problem was extended to irregular 
planar domains and solved by means of finite elements. 
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Weickert and Schnorr m considered a spatio-temporal model by extending 
the domain to x [0,T]. It additionally favours temporal regularity of the 
solution by including first derivatives with respect to time. Such models are of 
particular interest whenever trajectories are to be computed from the optical 
flow field. A unifying framework including several spatial as well as temporal 
regularisers was proposed in [JT]. For the purpose of evaluation and flow field 
visualisation a framework was created by Baker et al. [6]. 

Recently, generalisations to non-Euclidean domains have gained increasing 
attention. In m and m optical flow was considered in a spherical setting. 
Lefevre and Baillet [37] adapted the Horn-Schunck functional to surfaces embed¬ 
ded in Following Schnorr [35], they proved well-posedness of their formula¬ 
tion and employed a finite element method for solving the discrete problem on a 
triangle mesh. With an application to cell motion analysis, Kirisits et al. [221121I 
recently considered optical flow on evolving surfaces with boundary. They gener¬ 
alised the spatio-temporal model in [32] to a non-Euclidean and dynamic setting. 
Eventually, the problem was tackled numerically by solving the corresponding 
Euler-Lagrange equations in the coordinate domain. Similarly, Bauer et al. [7] 
studied optical flow on time-varying domains, with and without spatial bound¬ 
ary. They proposed a treatment on surfaces parametrised by product manifolds, 
constructed an appropriate Riemannian metric, and proved well-posedness of 
their formulation. 

In Kirisits et al. [23|, the authors considered various decomposition mod¬ 
els for optical flow on the 2-sphere. The proposed functionals were solved by 
means of projection to a finite-dimensional space spanned by vector spherical 
harmonics. Concerning projection methods, Schuster and Weickert m solved 
the optical flow problem in solely based on regularisation by discretisation. 

Regarding sphere-like surfaces and spherical harmonics expansion of closed 
surfaces we refer to [32] and the references therein. 

Einally, let us mention m [211131131], where optical flow was employed for the 
analysis of cell motion in microscopy data. In particular, in Schmid et al. [33] the 
embryo of a zebraflsh was modelled as a round sphere and motion of endodermal 
cells computed in map projections. 

The remainder of this article is structured as follows. In Sec. we formally 
introduce evolving sphere-like surfaces, recall the definition of vectorial Sobolev 
spaces on manifolds, and discuss both scalar and vector spherical harmonics 
on the 2-sphere. Section is dedicated to optical flow on evolving surfaces 
and our variational formulation. In Sec. |3] we discuss the numerical solution. 
In particular, we propose to solve the resulting energy in a finite-dimensional 
subspace and rewrite the optimality conditions to be defined solely on the 2- 
sphere. Moreover, we show how to fit a sphere-like surface to the labelled cells 
in the microscopy data. Einally, in Sec.[^ we solve for the optical flow field and 
visualise the results. The appendix contains deferred material. 

2 Notation and Background 

2.1 Sphere-Like Surfaces 

Let 

52 = {:r G RV ||rc|| = 1} 
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be the 2-sphere embedded in the 3-dimensional Euclidean space. The norm of 
n = {2,3}, is denoted by \\x\\ = By 

X : C ^ (3) 

we denote a smooth (local) parametrisation of 5^ mapping coordinates ^ = 
G O to points X = G S^. 

Furthermore, let I := [0, T] C M denote a time interval and let Ad = {Ait}tei 
be a family of closed smooth 2-manifolds Mt C Each t G /, is assumed 
to be regular and oriented by the outward unit normal field N(t,x) G 
X G Mt- We assume that M (locally) admits a smooth parametrisation of the 
form 


H. G M (4) 

and call M an evolving sphere-like surfaee. 

We denote by / : At ^ M a smooth function on the moving surface. Its 
coordinate representation / : / x ^ M and its corresponding spherical repre¬ 
sentation / : / X 5^ ^ M are given by 

0 = x(C)) = f{t, y{t, C)). (5) 

As a notational convention we indicate functions living on 5^ with a tilde and 
functions on M with a hat, respectively. Their corresponding coordinate version 
is treated without special indication. 

For convenience, we define smooth extensions of / and / to \ {0} by 

A*’» / (‘’ h) ^ G (*' h) h) ' 

respectively. Note that, while / is constant in the direction of the surface 
normal of 5^, the extension / in general is not. We point at Fig. [^illustrating 
the setting. 

Similarly, for vector-valued functions u : I X S‘^ ^ and u : M ^ R^ 
the extensions to R^ \ {0} are defined component-wise and for all times tel. 
They are denoted by u and u, respectively. As a notational convention, boldface 
letters are used to denote vector fields. Moreover, we distinguish between lower 
and upper case boldface letters. The former identify tangent vector fields and 
their extensions to R^ \ { 0 } whereas the latter indicate general vector fields in 
mA 

For a differentiable function / : / x ^2 ^ M, we write dif as an abbreviation 
for the partial derivative of / with respect to That is, Vm 2 / = (9i/, ^ 2 /)^, 
where Vm 2 is the gradient of M^. 

The tangent plane at a point y(t,<^) G Mt is denoted by Ty(^t,^)Mt and 
the tangent bundle by TMt = {{y(^:0} ^ ^ ^ orthogonal 

projector onto the tangent plane T^Mt at x G Aft, t G /, is given by 

P^(t,x) =Id-N(t,x)N(t,x)“^ 

In particular if Mt = 5 ^, that is p in 0 is identically one for all t e the 
outward unit normal and the orthogonal projector are given by N and P 52 , 
respectively. 
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Figure 3: Schematic illustration of a cut through the surfaces and Aft inter¬ 
secting the origin. In addition, a radial line along which the extensions / and / 
are constant is shown. The surface normals are shown in grey. 


In what follows, we define spatial differential operators. As they are identical 
to those on static surfaces we consider time t G / arbitrary but fixed. Then, the 
surface gradient of /, as given in ([^, is defined by 

VA</:=PAiVK3/eM3, (7) 


where Vrs denotes the usual gradient of the embedding space. Let us stress 
that it is independent of the chosen extension, see e.g. dSl p. 389]. 

We emphasise that, in particular, if Aft = for alH G / it follows that 


Vrs/ = P^sVrs/ + (Id - 


The last term of the sum on the right hand side is the normal derivative of /, 
which according to the definition of the extension in (§1 vanishes. Thus, 

V 52 / = P 52 Vrs/ = Vrs/. ( 8 ) 

For convenience let us observe that, by taking 9^/ in ([^, we arrive at 

= VM3/(t,x(0) • 9tx(0 = V 52 /(t,x( 0 ) • 9tx(0 (9) 

due to the chain rule and the projection onto the tangent plane Tx(^)iS^. 

Analogously to the surface gradient we define the spherical Laplace-Beltrami 
of / : / X 52 ^ M as 

A52/=-Ar3/, (10) 

where Ar 3 is the standard Laplacian of 
The set 

d2y{t, C} c (11) 

where y is the parametrisation defined in Q , forms a basis of the tangent space 
,t)Aft at y(t,^). Its elements form the gradient matrix Dy, which is derived 
as follows. 
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Let p be the extension of p : / x 5^ ^ (0, oo) according to Then, y 
from Q can be rewritten as 


y(^,0 = p(Lx(0)x(0- 


By the chain rule, 


diy{t,0 = (VM3p(i,x(e) ■ + p{t,x{ 0 )dix{ 0 . 

Using (|^ and the fact that p equals p on gives 

diY = (V52/5 ■ 9 jx)x + pdiX, 

where we have omitted the arguments (t, and (^) for better readability. When¬ 
ever convenient and no confusion will arise we will continue to do so. 

By applying ^ backwards and the fact that p(t, x(^)) = p(^) we have shown 

Dy = {diy d2y) 

= ((dip)x (d 2 p)x) + pDx e 


where Dx = (9ix, d 2 x) is the gradient matrix associated with x. 

As a consequence, we can uniquely represent a tangent vector u 
as u = ^^diy, where u = {u^,v?‘p € is its coordinate representation, 

see e.g. [26[ Prop. 3.15]. We call the components of u. 

In the sequel we will use Einstein summation convention. We sum over every 
index letter that appears exactly twice in an expression, once as a sub- and once 
as a superscript. For instance, we write u = u'^diy for the sake of brevity. 

We underline that the coordinate basis is not orthogonal in general. We 
will, however, require an orthonormal frame {ei(t,^),e 2 (t,^)} of the tangent 
space from Sec. 


2.2 


onwards. In the coordinate basis it reads 


= aidjY, 


(13) 


where : / x ^ M, i, j = { 1 , 2 }, are functions obtained from the Gram- 
Schmidt process. 

Combining <§ and with the expressions derived for Dx and Dy we can 
conveniently state that 


= DxJV and 


= Dy^VMf^ 


(14) 


Let us derive the following useful generalisation of For a tangent vector 
V = v'^diX G X G 5^, the directional derivative of / along v at x is 


Vfis/- V = V 52 / -u'^jX = (Dx'^Vs^f) • V = Vr2/ • V = v^dif, 


(15) 


where the third equality follows from the first equation in (14). Analogously, 
for V = v^diY ^ with x G Alt and t G I, one can derive 

Vx/-v = 'c*5i/. (16) 

As soon as we have established the relation between v and v it will conveniently 
allow us to switch between (15) and (16). 
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Moreover, the coordinate representation of the surface gradient is derived 
as follows. Let us start out with the first equation in (14). By writing V 52 / in 
the coordinate basis, that is V 52 / = I^xu for some u, we obtain from (14) 

Vm 2/ = L)x^I)xu. 

Multiplying with Dx.)~^ from the left yields 

Dx.)~^'VT^ 2 f = u. 


Thus, 

V 52 / = Dxu = Dx(Dx'^L>x)-Wr 2 /. (17) 

Furthermore, let t G / be fixed and let /(t, •) : Mt The surface integral 

of / is 

[ fdMt=[fJyd^^ (18) 

JMt dn 

where Jy is the Jacobian of y. According to Theorem 3 in m p. 88 ], it is given 
by 

{Jyf = det(Dy^Dy) 

and by using ( | 12 | ) yields 

{Jyf = y {{dipfd2X ■ d2X + {d2pfdix ■ dix 

+y (9ix • c1ix)(92X • d2x) - 2dipd2p{dix ■ d2x) - p^{dix ■ 32x)^). 

Note that x • x = 1 and thus, terms of the form 9^x • x vanish. By the differen¬ 
tiability of X, ones has 

di{x • x) = 0. 

Therefore, S^x-x = 0, meaning that tangential and normal vectors are orthogo¬ 
nal. We emphasise that Dy^ Dy is commonly referred to as Riemannian metric. 
It is positive definite and thus, (Ty(t,> 0 for all {t^C) G / x 11. 

The parametrisations x and y defined in ^ and respectively, suggest 
the straightforward construction of a smooth map 0 (t, •) : 5^ ^ H is given 
by the composition (y o x“^)(t, •), that is 


(/)(t, x) : X p{t, x)x. 


The differential D(j){t^x) : x)^t of 0 is a linear map and is given 

by 


D(j){t^ x) = p(t, x)Id + xVs^p{t^ x) G 


])3x3 


( 20 ) 


It follows from a direct calculation akin to the derivation of Dy in ( 12 ). 

Let us exhibit the action of x), for x = x(^) and t G /, onto a tangent 

vector V = t’^Q^x G TxS‘^. We have 


D^{t, x)(v) = p(t, x)v + x{\/s^p{t, x) • v) 

^p(t,x)v + xv^dip(^) 

= p{t, x)v^diX + xv^dip{^) ( 21 ) 

= u® {p{t, x)dix + xdip{e)) 
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In other words, the components are preserved whenever a tangent 


vector is mapped from S‘^ to Mt via the differential (20). 

As a matter of fact, given a tangent vector field v = v^diyi on 5^, the 
differential gives rise to a unique tangent vector field v = v^diy on 
see [ini Chapter 8]. Whenever we use v and v in the sequel we refer to their 
unique identification via the differential (20) and call v the pushforward of v. 


At this point, the reader might find it helpful to have a look at Fig. 

With the above definitions at hand we are able to relate the surface inte¬ 
gral (18) to an integral on 5^ via a change of variables. The key is to compute a 
meaningful surface element as \det{D^)\ is the magnitude of the change of the 
volume element. The following lemma provides the required form. 

Lemma 1. Let x : [0,7r] x [0, 27r) be the standard parametrisation of , 

^ (sincossinsin^^, cos^^)^, 

and let f : Ad ^ M and p : / x 5^ ^ (0, oo) he as above. Then, for t ^ I, 

f fdMt = [ fp./WTM^T^ds\ 

JMt ds‘^ 

Proof Let us denote by ei(^) and e 2 {f) the orthogonal unit vectors on 
in direction of and respectively, which are obtained by normalising the 
coordinate basis {9ix(^), 92x(^)}. That is. 


ei(0 = and §2(0 = 


^ 2 x (0 

ii52x(oir 


( 22 ) 


Moreover, a straightforward calculation gives 


Dx^Dx = 

and thus, the surface gradient of p in 


(1 0 \ 

\^0 sin^^^y 

spherical coordinates 0 is given by 


Vs2p(i,x(e) = dipi^dixiO + ■ I d2p{0d2x{0 

sm t, 

^dipiO ei(^) + -r^d2p{0 62 ( 0 , 

sin^^ 

where we have replaced the coordinate basis with the orthonormal basis. 

Using Dx^Dx in ( [l^ , the Jacobian Jy can be written as 

(^y)' = P {{dip? sin2 0 + {d2P? + y sin2 0) 

= P {{dip? + -rY^{d2P? + P) sin^ 
sm t, 

= P(l|V5^p||^+y)sin^0- 

Here, we have omitted the argument (t, x(^)) of Vs‘^P‘ Then, the integral turns 


10 










out to be 



^27r ^TT 

/ / fJy 

Jo Jo 

r27r pTT 

/ / fpV\\^s-p¥+P^ sine 

Jo Jo 

[ fpVW^W^dS^ 

Js^ 


where the last equation follows from (18) if Alt 
and 


Jx = y det{Dx.^ Dx.) = 


= 5^, the fact that sin^^ > 0, 
sin^^ 


□ 


The concepts introduced above, and further properties thereof, may be found 
in any standard differential geometry book. For instance, in [niiinii25i[26]. 


2.2 Vectorial Sobolev Spaces on Manifolds 

We briefly introduce the appropriate function spaces required for the variational 
optical flow formulation on Riemannian manifolds. Again, let us consider time 
t G / arbitrary but fixed. 

For a tangent vector field v on Aft we denote by the covariant 

derivative of v at x G Aft along the direction of a tangent vector u G 
We define it as the tangential part of the usual directional derivative of the 
extension v along u in the embedding space, that is, 

V(iv := Pa4Vr3v(u). (23) 


It is a linear operator Vv(x) 
is given by 


TxMt —> TxJAt and its Hilbert-Schmidt norm 


l|Vv(x)||i = y^llVe.vG) 


(24) 


i=l 


where {ei, 62} denotes the orthonormal basis of the tangent space T^^Aft, cf. (13). 
We stress that ( [^ is invariant with respect to the chosen parametrisation. 

For each t G /, we define the Sobolev space i^^(Aft, TAft) as the completion 
of the space of vector fields C^(Alt,TAft) with respect to the norm 




-j 

JMt 


lIVvIli dMt 


(25) 


where the surface integral is defined in (18). Let us emphasise that (25) is indeed 
a norm whenever Alt is diffeomorphic to the 2-sphere. The reason is that, by 
virtue of the Hairy Ball Theorem, there is no covariantly constant tangent vector 
field but V = 0, see e.g. m p-125]. 

Alternatively, one can define Sobolev spaces of vector fields such that each 
component of a vector field originates from a scalar Sobolev space. See, for 
instance, Lefevre and Baillet m- On the 2-sphere, however, they are typically 
introduced by means of the spherical Laplace-Beltrami operator, see e.g. m 


Chapter 6.2] and Sec. 2.3 for the scalar counterpart. For a thorough treatment 
of Sobolev spaces on Riemannian manifolds we refer to the books [a EH]. 
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2.3 Spherical Harmonics 

We denote by Harm^ the space of homogeneous harmonic polynomials of degree 
n G No with their domain restricted to S^. Its dimension is 


dim(Harm^) = 2n + 1. 


An element W ^ Harm^, n G No, is called a (scalar) spherical harmonic. It 
is an infinitely often differentiable eigenfunction of the Laplace-Beltrami opera¬ 
tor A52, defined in (10), with corresponding eigenvalue 

An = n{n + 1). 


We refer to Theorem 5.6 and Lemma 5.8 in [301 Sec. 5.1] for detailed proofs of 
the previous statements. 

The set 

: n e No,i = 1, ... ,2n + l| (26) 

is a complete orthonormal system of with respect to the inner product 

(•, •)^2(52) on 5^. In further consequence, for a function / G 1/^(5^), we have 
the Fourier series representation 


oo 2n+l 

n=0 j=l 

Again, we refer to isni Sec. 5.1] for the proofs, in particular to Theorem 5.25. 
In the present article we employ fully normalised spherical harmonics. For the 
explicit construction see [30l Sec. 5.2]. 

Moreover, the norm of L‘^{S‘^) is readily stated in terms of the coefficients 
in the above expansion via Parseval’s identity 


Il/Ili2(52) — ^(/: ^ni)i 2 ( 52 ). 

n,j 

For an arbitrary real number s G M, we define the Sobolev space as 

the completion of all (7^(5^) functions with respect to the norm 


11/11^3(5.) := ii(A5. + = E(^" + ryf,ynj)his^y 

n,j 


We stress that, by (10), A 52 / = — Aj^s/ and we have > 0 for all n G No 
yielding a sound definition. Accordingly, for 5 G M, we define the seminorm 
of order s by 


•— 11^5^2 f\\‘L^(S^) — ^n(/5 ^ni)i2(52y (27) 

n,j 

Now that the space L‘^{S‘^) is endowed with a basis, we can proceed to define 
an orthonormal system for square integrable tangent vector fields on the sphere. 
This will immediately allow us to treat vector-valued problems consistently. 
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Let Yn G Harm^ be a scalar spherical harmonic of degree n G No- Any 
vector field y : 5^ ^ that can be written in the form y = , where 

:= YnN, 
yi^) := 

yi^) := V52f„ X N, 


is called a vector spherical harmonic of degree n and type i, cf. [121 Defini¬ 
tion 5.2]. Recall that N is the outward unit normal to S^. 

By definition, yl^^ is a normal field whereas yl^^ and y^f^ are tangent vector 
fields. Consequently, the latter are called tangential vector spherical harmonics. 
Note that, by means of the Hairy-Ball Theorem, no tangential vector spherical 
harmonics of degree zero exist. 

In further consequence, let us denote by 1/^(5^,T5^) the space of square 
integrable tangent vector fields on equipped with the inner product 


(u,v)i^2(52^^52) = / U- 

Js^ 


vdS\ 


Here, dS^ denotes the usual spherical surface measure, see also Lemma [T] 

Since (26) is an orthonormal set for 1/^(5^), the set 

{yS :neN,j = l,...,2n+l,i = 2,3}, (28) 

is an orthonormal system for 1/^(5^, T5^), where we have defined 


yS = 

yi? = X N, 


(29) 


for orthonormalisation purpose, see [121 Sec. 5.2]. Thus, every vector field v G 
1/^(5^, TtS^) can be written uniquely as 


3 oo 2n+l 

^ = EEE (v,yi*j)L2(s2,rs2)yi*]. 

i=2 n=l j=l 

We refer to the books [121 EO] for further details on the subject. Table 
contains a summary of notation used in the sequel. 


3 Optical Flow on Evolving Surfaces 

3.1 Generalised Optical Flow 

Optical flow models are typically based on the assumption of constant brightness. 
Given a sequence of (planar) images 

y I / X (G M X — y M 

such that / G C^(/ X fl), it assumes that the intensity /(t, 7(t, ^)) stays constant 
over time when moving along a trajectory 7(',0 : / ^ D starting at ^ G 11. In 
other words, in the planar setting, we have 

^/(f, lit, O = dtf + Vr2/ ■ dti = 0, 
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n 

coordinate domain 

I 

time interval 

52 

2-sphere 

M 

family of sphere-like surfaces Mt 


tangent plane at x G 5^ 

TyMt 

tangent plane dX y ^ Mt 

N,N 

outward unit normals to and M 

X, y 

parametrisations of S‘^ and M 

Dx, Dy 

gradient matrix of x and y 

{ 5 iX, d2x} 

basis for TS‘^ 

{ 5 iy, d2y} 

basis for TM 

{ei, 62} 

ortho normal basis for TMt 

V 

surface velocity of M 

(i>,D(j> 

smooth map from to M and its differential 

f’J 

scalar function on S^, M, and their coordinate version 

Vx/ 

surface gradient on and Mt 

V, 'Cf, V 

tangent vector fields on M, and their coordinate version 

VfiV 

covariant derivative of v along direction u on Aft 

Ynj 

scalar spherical harmonic of degree n and order j 

J nj 

vector spherical harmonic of degree n, order j, and type i 

^ (i) 

Tnj 

pushforward of y^j via the differential 


Table 1: Summary of notation used throughout the paper. 


which is termed optical flow equation and must hold for all ^ G 17 and all t G /. 
For the sake of consistency, we denote by dt the partial and by d/dt the total 
derivative with respect to time. 

It is possible to generalise the idea to a non-Euclidean setting where the 
image lives on a, potentially moving, manifold. To this end, let us be given an 
evolving surface 

M:=[j ({t} X Mt) C 

tel 

specified by a parametrisation y : / x O ^ as in 0 together with a function 
/, its domain being M. For a time t G /, 

/(t, •) : Mt ^ 

is then an image on the surface. Adapting the above idea of constant brightness 
to the new setting requires that, along a smooth trajectory 7(', x) : t ^{t,x) G 
Mt that starts at x G Mq and always stays on the surface, we must have 

= f{0,x). (30) 

However, in order to proceed as above one needs to define a meaningful deriva¬ 
tive with respect to time. 

One possibility, which is pursued in 1221211 , is to consider derivatives along 
trajectories following the moving surface. Let y be as above and let dty = V be 
the surface velocity, its domain being ^ emphasise that 
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V is in general not tangent to Att, t G /, and hence in our notation denoted by 
a boldface capital letter. Then, 


dYf{to,xo) := 


(31) 


t=to 


is the time derivative of / at xq = y(to, 0 along the parametrisation y(', 0- 
a consequence, one can deduce that 


dYf = d^f + VMf-y 


holds, where df^f{to,xo) is the time derivative of / in normal direction. It is 
defined analogously to ( [M| albeit following a trajectory through xq G Aito 
for which xq) is orthogonal to Tx^Aito- 

From that one can immediately formulate the above idea of constant bright¬ 
ness (30) along 7. To this end, we define by M := 9^7 the velocity of a point 


moving along the trajectory 7 and demand that 


df/ = dN/ + VA,/-M = 0 


(32) 


must hold. Equation (32) is a generalised optieal flow equation. In Fig. we 
sketch the various trajectories through the evolving surface and their corre¬ 
sponding velocities. 

Since we are, however, interested in a coordinate representation of 7, we 
define a family of trajectories (3 : I x Q ^ Cl such that 


7(i,y(o,e) = y{t,l3{t,0) 

holds for alH G / and all ^ G In other words, we want the composition of [3 
with y, and 7 to coincide. By taking the total derivative d/dt on both sides of 
the above equation we get 


dtl = dty + dtp^diy. 


Let us denote v := dtjd^diy and recall that dty = V is the surface velocity. 
The above relation states that the total velocity M = 9^7 along a level line of 
constant intensity is 

M = V + V (33) 

and V is a tangential velocity relative to the prescribed surface velocity V. 


Solving the generalised optical flow equation (|32|), however, is inconvenient 
as 


and, in further consequence, df^ is unknown or hard to estimate. Never¬ 


theless, one can relate (32) and (33), as shown in El Lemma 2], and arrive at 


the parametrised optieal flow equation 


7/ +Va^/-v = 0. 


(34) 


Solving for the optical flow then means finding a (time-varying) vector field v 
that is tangent to the surface at all times and satisfies the above equation at 
every point x G At on the moving surface. 

Let us conclude this subsection with a remark that, in general, there exist 
infinitely many parametrisations y for a given evolving surface. The actual 
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Figure 4: Illustration of trajectories through the evolving surface. Their corre¬ 
sponding velocities are shown in grey. 


surface velocity however might be unknown or cannot be estimated from the 
data, as it is the case in this work. As a remedy we impose a surface velocity by 
choosing a natural surface parametrisation y. For a further discussion of this 
matter we refer to [7]. 

Moreover, we again stress that the sought tangent vector field v depends on 
the chosen y, or equivalently on V, and should be interpreted with care. The 
actual trajectories 7 though can be reconstructed by finding the integral curves 
of (33). For this precise approach we point the reader to m 


3.2 Variational Formulation 

The optical flow equation ( [M| ) derived above is underdeter mined and, in gen¬ 
eral, a unique solution is not ensured. A common technique to deal with non¬ 
uniqueness is Tikhonov regularisation, where one finds a minimiser of 

IMf/ + V^,/•v||7^)+a7^(v). 

Here, 7^(v) is a regularisation functional and a > 0 a regularisation parameter, 
balancing the two terms. The first term is typically referred to as data term 
whereas the second is called smoothness term. The latter enforces uniqueness 
and incorporates prior knowledge about favoured solutions. 

A common choice for 7^(v) is the squared Sobolev seminorm, involving 
first derivatives with respect to space and time. It favours spatial as well as tem¬ 
poral regularity and is of particular interest when trying to estimate trajectories 
of objects, albeit computationally more demanding. See, for example [snug 
and [7]. 

Alternatively, one can omit temporal regularisation leading to a regulariser 
of the form 

■^(v) = j\\Ht,-)\\‘H^Mt,TMt) dt, 

which is equivalent to solving for each time instant separately. It resembles the 
original formulation in m and its extension to 2 -manifolds m In the present 
article we follow this approach and attempt in finding the unique minimiser 
V G TAit of the energy 

fa(v) ;= \\dY f + ^Mf (35) 
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Figure 5; Commutative diagram relating spaces Q,, S^, and Mt, and tangent 


vector fields. We highlight that jp is the coordinate representation, see Sec. 2.1 


of a particular tangential vector spherical harmonic and is its uniquely 
identified tangent vector field on Alt. 


for each time instant t G I separately. Superimposing temporal regularisation 
however is straightforward, see m Sec. 2.2.2], but not considered here. 

4 Numerical Solution 


4.1 Finite-dimensional Projection 

For the subsequent discussion we let t G / be arbitrary but fixed and assume 
to be given a parametrisation y(t, •) : Q Alt as defined in We defer the 
question of how to find it to Sec. |4.3[ Moreover, for notational convenience, we 


relabel the set of tangential vector spherical harmonics (28) using a single index 


letter p G N. For instance, for the expansion of a tangent vector field on 5^ we 
simply write u = Upfp^ where G M are the coefficients. 

We intend to approximate the solution of the problem 

min £a (v) 


in a finite-dimensional subspace U C i^^(Alt, TMt), where £a is defined in (35). 
We define this space as 

U = spanjyp : p G J^/}, 

where C N is a finite index set and fp is the pushforward of a particular 
vector spherical harmonic fp via the differential see (20). Figure gives 
a descriptive view of the relation between the introduced spaces and tangent 
vector fields. 

The sought vector field is then uniquely expanded as 


peJu 


(36) 


with Vp G M, p G Ju, being the coefficients. Minimisation of functional (35) 


results in a finite-dimensional optimisation problem over Plugging ansatz 
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(36) into (35) and writing out definition (25) of the Sobolev {Ait Ait) norm 
gives 

^a(v) = / i{dY f+^Vp{yMf-fpif+ a\\^VpVypf^ dMt- (37) 

J Mt ^ T ^ 


V^Ju 


peJu 


By applying the definition of the Hilbert-Schmidt norm (24), using linearity 
of the covariant derivative V{iV with respect to v, and the definition of the norm 
of we obtain the representation 


iiv vpYpWi = y]ii Y ^'pW.ypH^ 

peJu *=i P^Ju 

2 

= Z!( Z! ■ Z 

i=i peJu q^Ju 

2 

= Z Z '^p^?(W.yp ■ Vg.y,) 

i=l p,qeJu 

for the regularisation term. 

The optimality conditions for the discrete minimisation problem (37) are 
obtained by taking dSajdvp = 0, for all p G Ju, and are given by 

r ^ 

•yp)(Vx/-y,) +aZ(Veiyp • Ve.y,)) dMt 

= - [ dfpVMf ■ fp) dMt, P e Ju- 

JMt 


q^Ju 


(38) 


In matrix form they read 

(^4 + aD)v = 6, 

where = (i;i,..., G is the vector of unknowns. The entries of the 

matrix A = {apq)pq are 


= [ {'^Mf ■yp){'^Mf-yq) dMt, 

JMt 


the entries of the matrix D = {dpq)pq associated with the regularisation term 
are given by 


dpq 


/Mt i=l 

and the entries of the vector b = {hp)p are 


JMt .-1 


bp = - [ dtf{VMf-yp)dMt. 
JMt 
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4.2 Rewriting the Optimality Conditions 


Even tough directly solving the derived optimality conditions (38) is perfectly 
legitimate, we take a different approach. The goal of this section is to rewrite 
the optimality conditions in terms of quantities defined on the 2 -sphere, thereby 
allowing a more general treatment. On the one hand, we want to deal with all 
surfaces M.t for all t G / in a unified manner and, on the other hand, we aim 


at evaluating (38) numerically on the (approximated) sphere without having to 
deal with multiple charts, see e.g. [9]. 

The following is a straight-forward generalisation of [231 Lemma 2 ]. 

Lemma 2. Consider time tel arbitrary but fixed. Let v = and 

V = v^diY tangent veetor fields on and A4t; respeetively, sueh that 

they are related via the differential (20). Then, the parametrised optieal flow 


equation (34) is equivalent to 


9t/ +V52 /-v = 0. 


Proof. According to the definitions (31) and ([^, we have 


dY fit,y{t,0) = ffit,yit,0) 

= dtfit,x{^)) 


and it remains to show the identity 

Va4/'V = V 52 /-V, 


where we have omitted the arguments (t, y(t, ^)) on the left and (t, x(^)) on the 
right hand side, respectively. It follows directly from the coordinate representa¬ 
tion of the directional derivatives (15) and ([T^. □ 


In order to give coordinate expressions for the terms in ( [38| ) arising from the 
regularisation term we locally choose an orthonormal frame {ei(t, ^), e 2 (t, ^)} 
of the tangent space, see (13). As a consequence, the sought tangent vector field 
V can be written as 

V = w^Bi (39) 

for some components . The reason for expressing the unknown in an 

orthonormal frame, rather than the coordinate frame, is to simplify matters 


with regard to the Hilbert-Schmidt norm (24) of the covariant derivative. 


However, the chosen Galerkin method expands the unknown v in terms of 
the pushfoward of vector fields which are defined on the 2 -sphere, cf. We 

necessarily need to establish the relation between the intended form (39) and 
the expression in terms of the coordinate frame. 


Lemma 3. Again, lett e I be arbitrary but fixed and let u = u^di^ be a tangent 
veetor field on . Then, for a tangent veetor field v = w^Bi on A4t, we have 
V = D^(vi) if and only if . 
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Proof. (<^) First, note that aj(a = 5ji. Expanding v gives 

V = w^Bi = w'^oP^djY = {a~^)\u^ oP^djY = u^djy = D^{u), 

where we have used ( [21] ), cf. also Fig. [^ 

(^) Suppose V = D(j){u). Let us take the inner product with on both 
sides. For the left hand side we have 


V • e. 




6 ji = w\ 


For the right hand side we first observe that, by inversion of the matrix a in (13), 
it holds that djy = {a~^)jB£. Then, 

D^(u) • ei = u^dey ■ ei 


= u^{a 

= u^{a~^y^6. 




• e,- 


= u^{a 

and we conclude that as required. 

With the above relation at hand we obtain the following form. 


□ 


Lemma 4. Let tel and let u = u'^diy and v = v^diy be two tangent veetor 
fields on Ait- Then, we have 

2 

Ve,U • Ve^v = '^DiUWiV^, 

i=i 

where 

DiU^ ;= a’;dk{{a~^)lu^) + i,j = {1,2}, 

and DiV^ are defined aeeordingly. 

denote the Christoffel symbols with regard to the ortho normal frame 
{^ 1 , 62 } and are defined as 

V^.e/c = T^j^Bj. 

We refer to [211 Lemma 3] for their derivation. 


(40) 


Proof. First let us show that, for u = w^Bj as in (39), it holds that 

V^.U = DjU^Bn. 


By the product rule for the covariant derivative (23), 


W^.W^Bj = BjV^.W^ P W^V^^Bj. 


(41) 


Consider the first term of the sum and let be represented in the coordinate 

id 


basis as in (13). Then, 


Linearity of the lower argument of the covariant derivative with respect to 
C‘^(A4t) functions, cf. (16), yields 
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and by realising that is just the directional derivative (16) along dkY 

we obtain 

= a’ldkw^. 


Moreover, in the second term of the sum in (41) we use definition (40). Thus, 
by summing up all terms in ([4T| we obtain 




Applying the previous lemma gives coefficients Diu^ and DiV^ in the intended 
form. Finally, it remains to observe that 


V^.u • V^.v = DiU^Bj • DiV^Bj 
2 

= Y,DiuWiV^, 


since by definition • Bj = 6ij . 


□ 


Finally, by combining Lemmas m and we are able to express the opti¬ 
mality conditions (38) in terms of integrals on the 2-sphere. Thus, we arrive at 
the optimality conditions 


r ^ 

E / ((^5=/ • yp) ■ y.) + “ E DiviDivij p^\\vs^p\\^ + p^ds^ 

= - [ dtfiVs^f ■ fp) p^\\Vs^p\\^ + p^dS^, peJu, 
Js‘^ 

^ n (42) 

where p\/\\Vs‘^P\\^ + arises from the Jacobian (19), see also Lemma W 

The entries of the matrices A, D and of the vector 6 , respectively, are then 
given by 


apg = [ (V 52 / ■ fpKVs^f ■ fg) p^/W^WT^dS^ (43) 

Js^ 

f ^ 

dpg= Y. DiviDivi p^\\Vs-p\\^ + p^dS\ (44) 

bp = - [ dtfiVs^f ■ fp) p^\\Vs^p\\^ + p^dS^. (45) 

Js^ 


4.3 Surface Parametrisation 

In order to actually compute the above optimality conditions it remains to 
determine the radius p : I x S‘^ ^ ( 0 ,oo) in the presumed parametrisation (|^. 
Again, we continue the discussion for one particular but fixed time tel and 
drop the argument whenever convenient. 

Estimating p(t, •) : 5^ ^ (0, oc) is closely related to surface interpolation 
from scattered data. Given noisy data p^ and a parameter /3 > 0 , it amounts 
to finding the unique minimiser of the functional 

^siP) •= \\P - P^\\\‘2{S‘^) + (46) 
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where s > 0 is a sufficiently large real number, cf. definition (27). The first term 


penalises deviation from the observed data whereas the second term enforces 
spatial regularity of the solution. 

In practice, however, > 0 evaluations {p^{xi) : Xi G are given 

at pairwise distinct points on the 2-sphere. In our particular application these 
correspond to taking the norm in of pairwise distinct sampling points lying 
on the sphere-like surface A4t' 

p\xi) = ||a:;i||, {0}, i = 1,.(47) 

where Xi = Xi/||xi|| is the radial projection onto 5^. We again point the reader 
to Fig. 

Furthermore, before turning to the numerical solution of ( [4^ , let us briefly 
discuss the regularity requirements. In [7], the authors demand twice continuous 
differentiability for both the manifold M-t and the map y(t, •) to obtain well- 
posedness of the optical flow problem. By definition of the parametrisation (|^ 
we require that p(t, •) G (7^(5^). As a consequence of Theorem 2.7 in [151 
Chapter 2.6] regarding Sobolev embeddings, the space for s > 3 is the 

appropriate choice, i.e. C (7^(5^). 

Numerically, we approximate the solution of the problem 

min 

peHs(s^) 

by considering a finite-dimensional subspace Q C and point evalua¬ 

tions (47). In contrast to above, the space 




span{Tp :p G Jq}, 

where Jq C Nq again is an index set, is spanned by scalar spherical harmonics. 
The sought function is expanded as 

P — Pp^p^ 

peJQ 


where the unknowns are the coefficients pp G M, for p ^ Jq. Plugging into (46), 
applying definition (27), and taking dTjdpp^ for allp G Jq^ gives the optimality 
conditions 


N 


N 


qEJq 


i=l 




(48) 


Denoting by p = (pi,..., p| j^j)^ G the vector of unknown coefficients. 


the equations (48) can be written in matrix-vector form as 

(L -h PM)q = c. 

The entries of the matrix L = [Ipq) 


pqjpq 
N 


are 


ipq 




i=l 


the matrix M = diag(Af,..., is a diagonal matrix, and 


N 


Cp — ||l^(x^). 


i=l 


22 







4.4 Numerical Approximation 


Let us finally discuss the numerical solution of the optimality conditions (42). 
In particular, one needs to (approximately) evaluate the integrals (43), (44), 
and (45). Even though integrals on the 2-sphere can be computed exactly 
and quadrature rules exist up to a certain degree, see e.g. mm, we instead 
prefer to use a triangulation together with an appropriate quadrature. The 
reason is that numerical quadrature on the sphere would have to be of rather 
high degree to reproduce small details and features of the data, contrary to 
the chosen quadrature, which can easily be refined up to the desired precision. 
Finally let us mention that, for a more accurate evaluation of the integrals, one 
can introduce an intermediate (radial) map from the polyhedron to geodesic 
triangles. See e.g. [H Sec. 7.2], 

We use a polyhedral approximation S'^ = (V, T) of the 2-sphere S^. It is 
defined by a set V = {ui,..., C of vertices and a set T = {Ti, ..., T^} C 

V X V X V of triangular faces. Each triangle is most easily parametrised using 
barycentric coordinates, see e.g. [8j Chapter 5]. We associate with each triangle 
Ti e T di tuple identifying the corresponding vertices 

which are arranged in clockwise order. The parametrisation then reads 


MO = vh -ViO + OM -ViO 


with 

= {^ G G [0,1] and G [0,1 - ^^]}, 

which is referred to as the reference triangle. The gradient matrix of Ti is then 
simply 

Dxi = [diXi d2Xi) = (vi3 - Vii Vi2 - Ui J . 

The surface normal is constant on Ti and is denoted by N^. 

We approximate all functions on by corresponding functions on the poly¬ 
hedron 5^. A continuous function / : 5^ ^ M is replaced by its piecewise 
polynomial interpolation fh'.S‘^^RonS‘^. We define it as 


Nh 


fhi-) — 


(49) 


i=i 


Here, {(fj} are = 6 quadratic shape functions forming a nodal basis together 

with nodal po ints {vj} C 5^ and / is the usual radially constant extension. 
In other words, fh is both a radial projection from the 2- 


2.1 


cf. ^ in Sec. _ 

sphere to the polyhedron S‘^ and to piecewise quadratic functions. Note that 
the shape functions are defined on the triangular faces Ti. Whenever a function 
/ has a dependence on time we simply compute its approximation separately 
for all times tel. We point the reader to Fig. for a figurative illustration. 

In further consequence, the fully normalised scalar spherical harmonics, 
which were introduced in (26), are substituted with their corresponding ap¬ 
proximations on S‘^. For Y G Harm^, n G No , we have 


Nh 

n{-) = 


j=i 


(50) 
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Ph{j^Vl)vi 


Figure 6: Illustration of a triangular face (filled gray) intersecting the sphere 
at the vertices (hollow circles). The six nodal points consist of the vertices of 
the triangle together the edge midpoints (filled black dots). The approximated 
sphere-like surface is shown by the hatched gray area. A radial line passing 
through the vertex Vi is shown. The hollow circle indicates the i nters ection with 
S‘^ at which f{vi) in is taken. / itself, as described in Sec. |5.2[ is assigned 


by taking the maximum image intensity along the drawn radial line between 
the two cross marks. 


We chose piecewise quadratic approximations for Y so that we can ade¬ 
quately apply V c 2 and obtain piecewise linear vector fields. Accordingly, we 
define approximations of the vector spherical harmonics, introduced in (29), as 
follows. 


Proposition 1. Let Y G Harm^; n G N. The pieeewise linear interpolations 
of the eorresponding tangential veetor spherieal harmonies on a triangular faee 
Ti eT are 




i=i 


(?)), 




-v-1/2 Nh 

I *1 j=i 


(51) 

(52) 


Their derivation is deferred to the appendix. 

Without loss of generality, let //^(O, •) and //z,(l, •) be the approximations of 
the data / at two subsequent frames. We define the derivative with respect to 
time by the forward difference 

dtfh{-) ■■= A(i,-) - A(o, •)• 
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Moreover, we replace the surface gradient Vs‘^ f of a function on with its 
counterpart V 52 on , which is computed according to 0 - The function p 
is obtained by solving (48) and, for numerical computations, is further replaced 
with its piecewise quadratic interpolation ph as in ( [4^ . Coefficients are 
computed by the Gram-Schmidt process at the nodal points. For numerical 


computations piecewise quadratic approximations, as defined in (49), are used. 


Finally, for the calculation of the integrals we employ the standard quadra¬ 
ture on triangulated spheres, see e.g. [3l|T6]. Let = (1/3,1/3)^ be the centroid 
of the reference triangle fl. Then, we approximate the spherical integral over a 
function / : 5^ ^ M on the 2 -sphere by 

ppm 

/ fdS^^ fh dSl 

JS- JSI ^ 


5 Experiments 

5.1 Microscopy Data 

The present data consist of volumetric time-lapse (4-dimensional) images of a 
live zebrafish embryo during the gastrula period. These videos were recorded 
approximately five to ten hours after fertilisation by means of confocal laser¬ 
scanning microscopy and feature endodermal cells expressing a green fluores¬ 
cence protein. As a consequence, these labelled cells are recorded without back¬ 
ground and allow for a separate treatment. We refer the reader to m for many 
illustrations and a detailed discussion of the zebrafish’s developmental process. 
Regarding the imaging techniques used during data acquisition we refer to [28] 
and for the treatment of the specimen we point the reader to m- 

The crucial feature of endodermal cells is the fact that they form a so-called 
monolayer during early morphogenesis, see m- Essentially, it means that the 
labelled cells do not sit on top of each other but float side by side forming an 
artificial sphere-shaped layer. It can be regarded as a surface and allows for the 
straightforward extraction of an image sequence. Clearly, this surface is subject 
to geometric approximations. For instance, in [231311 it is assumed an ideal 
sphere, whereas in [7] and m only a fraction of the data is considered and 
modelled as a moving manifold and a height field, respectively, both possessing 
a boundary. 

The recorded data features a cuboid region of approximately 860 x 860 x 
320 pm^ of the animal hemisphere. The spatial resolution is 512 x 512 x 44 voxels 
and the recorded image intensities are in the range {0,..., 255}. Our sequence 
contains 75 images with a temporal interval of 240 s. For the further discussion, 
we denote the data by 


G {0,..., 255}'^^^^^^^^^^^^"^. 

5.2 Preprocessing and Surface Data Acquisition 

Let us briefly discuss the preprocessing steps required to obtain an image se¬ 
quence together with the evolving surface. We limit our consideration to two 
consecutive frames and denote the respective volumetric data by /q and ff. 
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Figure 7: Frames no. 70 (left) and 71 (right) of the processed image sequence 
in a top view. The embryo’s body axis is oriented from bottom left to top right. 


For each frame, the approximate surface is found by minimising the func¬ 


tional (46) with approximate cell centres acting as sample points. They appear 
as local maxima in image intensity and are readily located by Gaussian filter¬ 
ing followed by plain thresholding. However, beforehand the points are centred 
around the origin by first fitting a sphere and subsequently subtracting the 
spherical centre. 

The triangle mesh 5^ is obtained by iterative refinement of an icosahedron 
that is inscribed in the 2-sphere, see e.g. [HI Chapter 1.3.3]. Every refinement 
step halves the edge lengths by connecting the edge midpoints and projecting 
them to the unit sphere. Consequentially, every triangular face is split into four 
smaller triangles and the total number of faces after k G No subdivisions is 20*4^. 
In our case, k = 7 refinements are required to resolve the data adequately. 

It remains to discuss the acquisition of the approximations //z,(0, •) and 
//i(l,') on the polyhedron. For a frame j G {0,1}, we define the value at a 
nodal point Vi G 5^ in (49) via the projection 


f{j,Vi):= max jHcph{j,Vi)vi), 
ce[i-£,i+£] 


where 5 > 0 is chosen sufficiently large, fj denotes the piecewise linear extension 
of fj to which is necessary for gridded data. The above projection within 
the narrow band 

[(1 - ^)p/i(i, Vi)vi, (1 + e)ph{j, Vi)vi] 

corrects for small deviations of the cells from the fitted surface. Again, we 
refer to Fig. for illustration. Finally, all intensities are scaled to the interval 
[0,1]. Figure^ shows two frames of the extracted image sequence defined on 
the sphere-like evolving surface. Figure [^depicts the same matter but in a top 
view. For better illustration we have added an artificial mesh. Its radius has 
been widened by one percent. 
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5.3 Visualisation of Results 

We employ the standard flow colour-coding [5] for the visualisation of the com¬ 
puted vector fields. Its purpose is to create a colour image by assigning every 
vector a colour from a pre-defined colour disk. The colour associated is deter¬ 
mined by a vector’s angle and its length. 

However, it was originally defined for planar vector fields and requires adap¬ 
tation to our particular purpose of tangent vector field visualisation. To this end, 
we follow the idea developed in m by first projecting each vector to the plane 
and then rescaling its length. Let us denote by : (xi, X 2 , ^ 3 )^ (xi, X 2 , 0)^ 
the orthogonal projector of onto the xi-X 2 -plane. For a tangent vector field 
V we apply the colour-coding to the planar vector field 



It is constructed so that the length of individual vectors is preserved. Subse¬ 
quently, the obtained colour image is mapped back onto the surface. Clearly, 
in the above construction, one has to distinguish the cases where X3 > 0 and 
xs < 0. Moreover, P^^g is required to be injective in either case. 

The radius R of the colour disk is chosen to be equal to the longest vector 
in the respective vector field we attempt to visualise. Table lists all values of 
R for the different figures in this section. In Fig. we show a colour-coded 
tangent vector field together with the colour disk. 

For simplicity reasons, for image functions as well as surfaces we plot their 
piecewise linear approximations. Moreover, the visualised vector fields are eval¬ 
uated at the centroids and result in piecewise constant colour-coded images. 


5.4 Results 


We performed several experiments on said zebrafish microscopy data. In order 
to obtain an approximation of the evolving surface, we minimised functional 1 


by solving the optimality conditions (48). As mentioned in Sec. 5.2 approximate 


cell centres serve as input. The parameter of the Sobolev space H^{S) was 
chosen as s = 3 + e, where e = 2.2204 • 10“^^ is the machine precision, cf. also 
the discussion regarding theoretical requirements in Sec. |4.3[ The regularisation 
parameter was set to /3 = 10 “^ and the finite-dimensional subspace was chosen 
as 

Q = spanjWi : n = 0,... ,30, j = 1 ,..., 2 n + l|. 

In the second step, we computed a minimiser of functional as outlined 


in Sec. 4.4 Here, the finite-dimensional subspace was chosen as 


U = span|y^*j : n = 1,.. .,50, j = 1,..., 2n + 1, i = 2, sj. 


The linear systems resulting from optimality conditions ( [4^ and ( [48| ) were 
solved by means of the General Minimal Residual Method (GMRES) using an 
Intel Xeon E5-1620 3.6 GHz workstation equipped with 128 GB RAM. Solutions 


to (42) and (48) converged within 1000 and 100 iterations, respectively, to a 
relative residual of 10“^. The overall runtime was dominated by the evaluation 


of the integrals (43), (44), and (45). In our Matlab implementation it amounts 
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Figure 8: Function ph obtained by minimising for frames 70 (left) and 71 
(right). Colour corresponds to the radius of the fitted surface. S‘^ is depicted 
in a top view. 


Figure 

R 


10 11 (a) 11 (b) 11 (c) 11 (d) 12 (a) 12 (b) 12 (c) 12 (d) 


5.18 9.86 5.18 3.64 2.52 9.86 5.18 3.64 2.52 


Table 2: Radii R of the colour disks used for colour-coded visualisation of 
tangent vector fields. 


to several hours. However, the resulting linear system can typically be solved 
within seconds. Both implementation and data are available on our websiteQ 

Figures and portray a minimising function of for frames 70 and 
71 of the image sequence. The resulting surface is depicted in Fig. and in 
Fig.0 Clearly, it reflects the geometry appropriately and contains the desired 
cell features, cf. also the unprocessed microscopy data in Fig. 

In a second step we solved for minimisers of Eq, for different values of the 
regularisation parameter a. Figure [T^ depicts the optical flow field for a = 10“^ 
The tangent vector field is visualised as discussed in Sec. 


5.3 


Note that in 

all figures the colour disk has been scaled for better illustration. In Fig. El 
we illustrate tangent vector fields by means of the colour-coding obtained for 
a = 10“^, a = 10“^, (a = 1, and a = 10. Finally, in Fig. 12 we show the same 
results but in a top view. 


6 Conclusion 

With the goal of efficient cell motion analysis we considered optical flow on 
evolving surfaces. As a prototypical example we restricted ourselves to surfaces 
parametrised from the round sphere and showed that 4D microscopy data of a 
living zebrafish embryo can be faithfully represented in this way. In contrast 
to previous works, where only a section of the embryo or a spherical approxi¬ 
mation was considered, our approach fully attributes the geometry and models 
the embryo as as closed surface of genus zero. The resulting energy functional 
was solved by means of a Galerkin method based on vector spherical harmonics. 
Moreover, the parametrisation of the moving sphere-like surface was obtained 

^ http://www.CSC.univie.ac.at 
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Figure 9: Illustration of the function ph for frames 70 (left column) and 71 (right 
column). The bottom row differs from the top row by a rotation of 180 degrees 
around the xs-axis. 


from the data by solving a surface interpolation problem. Scalar spherical har¬ 
monics expansion allows to easily meet the smoothness requirements of the 
surface. Finally, we conducted several experiments based on said microscopy 
data. Our results show that cell motion can be indicated reasonably well by the 
proposed approach. 
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Figure 10: Tangent vector field minimising Depicted is the colour-coded 
optical flow field computed between frames 70 and 71. The right image differs 
from the left by a rotation of 180 degrees around the xs-axis. 


Appendix 


It remains to give the calculations regarding the piecewise linear approximations 
of vector spherical harmonics on S‘^. Both equations in Prop, [^follow directly 
by expanding the definitions of the tangential vector spherical harmonics (29). 
For the fist identity, that is (51), we have 


Nh 


yf = = A- 1/2 

i=i 


The second identity, that is (52), follows by the fact that 

ol-rl lo O I tCt X a2Xi 

2 Ti = X 52Xi , Ni = T- - --r, 

lOlXj X a2Xj| 

and by application of the vector triple product rule, yielding 


Th 


(3) = \-V^Vg,Yh X N,: 


= A-i/Wc2 V 


S‘^^h X 


9iXi X 92Xi 
l^lXi X a2Xi 


- 1/2 


217;, 


• a2Xi)aiXj - {Vsf^Yh ■ diXi)d2Xi), 


where the last equality results from the definition of the interpolation (50) of 
Yh and the directional derivative on the triangular face, that is 


Nh 

'^sNh ■ dkXi = {vj)dk‘Pj. 


i=i 
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Figure 11: Visualisation of the optical flow field obtained for different values 
of a. The bottom row differs from the top view by a rotation of 180 degrees 
around the xs-axis. From left to right: a) a = 10“^, b) a = 10“^, c) a = 1, 
and d) (a = 10. 
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a = 10. 
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